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HOMOGENIZATION  AND  CONTROL  OF  LATTICE  STRUCTURES 


G.L.  Blankenship^ 


March  2,  1985 


Abstract:  Under  certain  natural  conditions  the  dynamics  of  large,  low  mass  lattice  structures 
with  a  regular  infrastructure  are  well,  approximated  by  the  dynamics  of  continua,  e.g.,  trusses  may 
be  modeled  by  beam  equations.  Using  a  technique  from  the  mathematics  of  asymptotic  analysis 
called  homogenization,  we  show  how  such  approximations  may  be  derived  in  a  systematic  way 
which  avoids  errors  made  using  “direct”  averaging  methods.  We  also  develop  a  model  for  the 
combined  problem  of  homogenization  and  control  of  vibrations  in  lattice  structures  and  present 
some  preliminary  analysis  of  this  problem. 
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1.  Homogenization  of  regular  structures 


It  is  now  generally  accepted  that  large,  low  mass  lattice  structures,  e.g.,  trusses,  are  natural 
for  space  applications.1  Their  large  size  and  repetitive  infrastructure  require  special  techniques  for 
structural  analysis  to  cope  with  the  large  number  of  degrees  of  freedom.  Approximations  of  such 
systems  by  continua  provide  a  simple  means  for  comparing  structural  characteristics  of  lattices 
with  different  configurations,  and  they  are  effective  in  representing  macroscopic  vibrational  modes 
and  structural  response  due  to  temperature  and  load  inputs.  Our  approach  to  the  construction  of 
such  models  is  based  on  a  technique  for  asymptotic  analysis  called  homogenization.  It  has  been 
widely  used  in  mathematical  physics  for  the  treatment  of  composite  systems  like  porous  media  for 
which  one  wishes  to  have  an  effective  approximating  system  with  parameters  which  are  constant 
across  the  structure.2  Before  developing  the  general  features  of  the  method  and  applying  it  to  the 
treatment  of  lattice  structures,  we  shall  make  a  few  remarks  on  other  work  on  continuum  models 
which  has  appeared  in  the  recent  structural  mechanics  literature. 

Noor,  et.  al.  (1978)  use  an  energy  method  to  derive  a  continuum  approximation  for  trusses 
with  triangular  cross  sections  in  which  the  modal  displacements  of  the  truss  are  related  to  a 
linearly  varying  displacement  field  for  an  equivalent  bar.  Plates  with  a  lattice  infrastructure  are 
also  treated.  In  (Dean  and  Tauber  1959)  and  (Renton  1969)  exact  analytical  expressions  for  the 
solutions  of  trusses  under  load  were  derived  using  finite  difference  calculus.  By  expressing  the 
difference  operators  in  terms  of  Taylor’s  series  Renton  (1970)  was  able  to  derive  continuum 
approximations  to  the  finite  difference  equations  resulting  in  expressions  for  equivalent  plate 
stiffnesses,  for  example.  In  a  recent  paper  Renton  (198-1)  used  this  approach  to  give  equivalent 
beam  properties  for  trusses,  which  complements  the  earlier  work  of  Noor,  Anderson  and  Greene 
(1978)  and  Nayfeh  and  Hefzy  (1978).  (See  also  (Anderson  1981).) . 

In  most  cases  a  continuum  model  is  associated  with  the  original  (lattice)  structure  by 
averaging  the  parameters  of  the  lattice  over  some  natural  volume  (e.g.,  of  a  cell  of  the  structure) 

^ee,  for  example,  (Aswani  1082),  (Juang  1084),  (Mikulus  1078),  (Renton  1070),  and  (Taylor  and  Balakrishnar.  1084), 
and  the  references  therein. 

2See,  for  example,  the  papers  of  Larsen  (1075  1976),  Keller  (1077),  and  the  reports  of  Babuska  (1975)  for  applica¬ 
tions  and  discussions  of  design  techniques. 


and  identifying  the  averaged  parameter  value  (mass  density,  stress  tensor,  etc.)  with  the 
corresponding  distributed  parameter  in  the  continuum  model.  A  specific  form  for  the  continuum 
model  is  postulated  at  the  outset  of  the  analysis;  e.g.,  a  truss  with  lattice  structure  will  be  approx¬ 
imated  by  a  beam,  with  the  beam  dynamical  representation  assumed  in  advance.  While  this 
approach  has  an  appealing  directness  and  simplicity,  it  has  some  problems. 

First,  it  is  very  easy  to  construct  an  example  in  which  the  approximate  model  obtained  by 
averaging  the  parameters  over  a  cell  is  not  a  correct  approximation  to  the  system  behavior.  This 

is  done  in  the  next  section.3  Second,  one  cannot  use  this  procedure  to  obtain  “corrections”  to  the 
approximation  based  on  higher  order  terms  in  an  expansion,  which  may  sometimes  be  done  in  an 
asymptotic  analysis.  These  terms  can  be  used  to  describe  the  microscopic  behavior  (e.g.,  local 
stresses)  in  the  structure.  Third,  the  averaging  method  (averaging  the  parameters  over  space) 
does  not  apply  in  a  straightforward  way  to  systems  with  a  random  structure,  since  the  appropri¬ 
ate  averaging  procedure  may  not  be  obvious.4  Fourth,  the  method  cannot  be  naturally  imbedded 
in  an  optimization  procedure;  and  controls  and  state  estimates  based  on  the  averaged  model  may 
not  be  accurate  reflections  of  controls  and  state  estimates  derived  in  the  course  of  a  unified  optim¬ 
ization  -  averaging  procedure.  In  particular,  the  method  does  not  provide  a  systematic  way  of 
estimating  the  degree  of  suboptimality  of  controls  and  state  estimates  computed  from  the  ideal¬ 
ized  model. 

In  this  work  we  use  a  totally  different  technique  called  homogenization  from  the  mathemati¬ 
cal  theory  of  asymptotic  analysis  to  approximate  the  dynamics  of  structures  with  a  repeating  cel¬ 
lular  structure.  Homogenization  produces  the  distributed  model  as  a  consequence  of  an  asymp¬ 
totic  analysis  carried  out  on  a  rescaled  version  of  the  physical  system  model. 

Unlike  the  averaging  method,  homogenization  can  be  used  in  combination  with  optimization 
procedures;  and  it  can  yield  systematic  estimates  for  the  degree  of  suboptimality  of  controls  and 
estimators  derived  from  idealized  models.  While  our  results  are  stated  in  terms  of  simple  st.ruc- 

3See  the  numerical  experiments  In  (Bourgat  1978). 

htomogenization  methods  do  apply  to  systems  with  a  randomly  heterogeneous  structure,  see  (Papanicolaou  and 
Varadhan  1982)  and  (Kunnemann  1983).  We  shall  not  treat  such  systems  in  this  paper. 
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tures,  they  demonstrate  the  feasibility  of  the  method;  and  they  suggest  its  potential  in  the 
analysis  of  structures  of  realistic  complexity. 

In  section  2  we  give  an  example  derived  from  (Bensoussan,  Lions,  and  Papanicolaou  1978) 
illustrating  some  of  the  subtleties  of  homogenization,  particularly  in  the  context  of  control  prob¬ 
lems.  In  section  3  we  derive  a  homogenized  representation  for  the  dynamics  of  a  lattice  structure 
undergoing  transverse  deflections.  We  show  that  the  behavior  of  the  lattice  is  well  approximated 
by  the  Timenshenko  beam  equation;  and  we  show  that  this  equation  arises  naturally  as  the  limit 
of  the  lattice  dynamics  when  the  density  of  the  lattice  structure  goes  to  infinity  in  a  well  defined 
way.  The  problem  of  vibration  control  of  a  lattice  is  posed  and  discussed  in  section  4.  In  section 
5  we  derive  a  diffusion  approximation  for  the  thermal  conductivity  of  a  one-dimensional  lattice 
structure.  This  property  is  useful  in  analyzing  new  materials  for  large  space  structures.  An 
operational  calculus  for  homogenization  is  sketched  in  the  Appendix. 

Acknowledgements:  We  are  grateful  to  Professor  George  Papanicolaou  for  bringing 
Kunnemann’s  paper  to  our  attention  and  to  Drs.  A.  Amos  and  R.  Liudberg  for  their  comments 
on  an  earlier  version  of  this  work. 


2.  A  one-dimensional  example 

From  (Bensoussan,  Lions,  and  Papanicolaou  1978)  we  have  the  following  example: 

d  ,  ,,  ,  duc(x)  ,  ,  ,  ,  .  ,  .  . 

- r-[a(z)  — -]  =  /(*).  x  e  (x  0,x  j)  (2.1) 

ax  ax 

UC{x0)  =  0  =«'(!,) 

where  a '(a-)  =  a  (x  /e),  and  a(y)  is  periodic  in  y  with  period  Y0,  a  (y)  >  a  >  0.  It  is  simple  to 
show  that 


|i  £  /  l«f(*)|2  + 


dv  c(x ) 
dx 


|  2  dx  <  c 


(2.2) 
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and  so,  u(  — +  u  weakly  in  the  Hilbert  space  H1.6  Moreover, 

ac  — ►  M  (a  )  =  f  a{y)dy  (2.3) 

o 

and  it  is  natural  to  suppose  that  u e  — ►  u  with  the  limit  defined  by 

-  \M(a)  -j-  u{x)]  =  /  (z),  x  E  (x0-*i)  (2.4) 

u  (x0)  =  u  (a:,) 

This  is  untrue  in  general  (Bensoussan,  Lions,  and  Papanicolaou  1978,  pp.  8-10).  The  correct 
limit  is  given  by 

~  «(*)]  =  /  (x  ),  X  e  (Xq.Xj)  (2.5) 

u  (x0)  —  «(ar,) 

with 

a  £  [M(  —  )]-1  (2.6) 

a 

In  general,  M(o)  >  a";  and  so,  the  error  in  identifying  the  limit,  (2.4)  versus  (2.5),  is  fundamen¬ 
tal. 

The  system  (2.4)  corresponds  to  averaging  the  parameter  a\x)  over  a  natural  cell;  a  pro¬ 
cedure  similar  to  that  used  in  the  past  to  define  continuum  models  for  lattice  structures.  As  (2.5) 
shows,  the  actual  averaging  process  can  be  more  subtle  than  one  might  expect,  even  for  simple 
problems. 

2.1.  Homogenization  of  the  example 

To  see  how  (2.5)  arises,  we  can  use  the  method  of  multiple  scales  which  applies  to  a  variety 
of  problems.  Suppose 


hlere  H1  =  {  U  6  L2(x0,x  J  :  |  |«  | \Hi  <  oo  } 
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(2.7) 


U  ((x  )  =  UC(x  )  =  U0(x  )  +  £  U  i(x  ,  )  +  ■  ■  ■ 

e  e  £ 

that  is,  we  suppose  that  depends  on  the  slow  scale  x  and  the  fast  scale  y  = 
adopt  an  ansatz  which  reflects  this  dependence.  Using  the  identity 


_d_ 

dx 


du 

dx 


1  du 
e  dy ' 


y 


x 

e 


then  (2.1)  may  be  rewritten  as 


.8  1  8  8  id., 

(  -f  —  — —  )  -  {  a  (y  )  (  — —  +  —  —  )  [  m0  -f  £  w,  +  •••]} 

ox  e  ay  ox  e  oy 


Simplifying  and  equating  coefficients  of  like  powers  of  e,  we  find  first  that 


J_  d_ 

e2  dy 


a{y)  -r-  «0]  =  o. 
dy 


The  assumptions  on  a(y)  imply 


«o {*  <y)  =  «<>(*) 


i.e.,  no  y-dependence.  The  coefficients  of  e  1  satisfy 


a  (y . 
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d  ,  .  .  d 

-xr-  l  a(y)  -7T-  «o. 

cJx  dy 


8_ 

dy 


[ « (y )  « !  ] } 

dy 


or 


_d_ 

dy 


,  ^  d 

a  (y)  -rr~ 

dy 


da 

dy 


8u0 

~d7 


If  we  look  for  u  j  in  the  form 


(x,y) 


-  x(y) 


du0 

dx 


+  V  i(x)  , 


then  the  corrector  x(y  )  must  satisfy 


d_ 

dy 


a  ( V ) 


~  xiv )  1 

dy 


da 

dy 


and  be  periodic.  That  is. 


x  /e;  and  we 


(2.8) 


=  /  (2.9) 


(2.10) 


(2.11) 


0  (2.12) 


(2.13) 


(2.14) 


(2.15) 
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a{y)  —r-  =  a  (y)  +  c 
ay 


(2.16) 


which  has  a  periodic  solution  (unique  up  to  an  additive  constant  in  y )  if  and  only  if 


~  I  1 1  +  -fr-  1  dy  =  0 
y o  o  a(y ) 


(2.17) 


which  implies 


c  =  -  [M(-)]-1  *  a 
a 


We  obtain  an  equation  for  n0{x)  from  the  solvability  condition  for  u2(x  ,y).  Equating  the 
coefficients  of  e°  in  the  expansion,  we  have 


d  d  d  d 

-  -7-  i  a(y)  —  u  2  -  —  a(y)  —  w ,  ] 

ay  dy  dy  dx 

d  .  d  .  d  d 

a7  '"(!,)af  *■'  “  aJ  « o  J  =  /  (* ) 


(2.19) 


This  has  a  solution  uQ(x  ,y ),  periodic  in  y  if  and  only  if 

Y  o 

{  -T7-  /  l  a  (y  )  +  4-  [  a  (y  )  x(y  )  ]  -  a(y)  4-  x(y)}  dy  }  .  d  t  (2.20) 

loo  vy  dy  dx2 

+  /(*■)  =  o 

where  we  have  used  (2.14).  The  integral  of  the  second  term  is  zero,  since  it  is  the  integral  of  the 
derivative  of  a  periodic  function  over  one  period.  Using  (2.16)  and  (2.18),  (2.20)  reduces  to 


d2u  o 

-  a  7  +  /  (»)  =  0 
(lx  ~ 


(2.21) 


(plus  the  boundary  conditions)  which  is  (2.5)  (2.6). 


2.2.  Control  and  homogenization  of  the  one  dimensional  system 

One  of  the  simplest  stochastic  control  problems  associated  with  the  preceding  system  is 
defined  by  the  Hamilton  -  Jacobi  -  Bellman  equation 


6 


,*  >  d°u< 

a  (— )  -r~r 

€  dx 2 


e  e  dx 


inf  [—  v 2  4-  g  (x  ,y)  v 
v  e  R  2  az 


-cm 


(2.22) 


G  O,  «{(z)  =  0  on  T  £  30 


where  O  is  an  open  interval  in  R,  and  each  function  a{y),  b(y),  and  g  (x  ,y)  is  periodic  in  y 
with  period  Y0.  We  assume  that  a  (y )  >  a  >  0  and  that  c  >  0,  and  that  the  controls  v  take 
values  in  R. 

This  Bellman  equation  corresponds  to  the  stochastic  control  problem 


u  e(x  )  =  inf  Je[v  (.)] 

t>  (•) 


J‘[v{.)}  =  E,  {  f  /(*',—  ,v)  [ 
0  6 


-  f  c  (ie,—  ,v)  it 


)dt  } 


(2.23) 


dx((t  )  =  cr(a:<.—  )  dw  (t )  +  —  b  (x( )  dt  +  G{xc,—  ,v)  dt 
e  e  e  e 

x((0)  —  x  GO,  t  >  0. 

with  cr{x  ,y  )  =  a  (y),  b(x,y)=b{y),  G  (x  ,y  ,v  )  =  g  (x  ,y  )v  ,  /  (x  ,y  ,v  )  =  ~  v2,  and 

c  (x  ,y  ,v  )  =  c  ,  a  constant  in  (2.22).  Each  function  in  (2.23)  is  assumed  to  be  periodic  in  y  with 
period  one.  We  are  interested  in  the  behavior  or  the  optimal  cost  and  control  law  for  (2.22)  in  the 
limit  as  e  — +  0.  The  stochastic  control  problem  (2.23)  was  treated  in  (Bensoussan,  Boccardo,  and 
Murat  1984);  the  analysis  here  uses  different  arguments  which  emphasize  the  computational 
aspects  of  the  system. 

Evaluating  the  infimum  in  (2.22),  we  have  the  nonlinear  system 


,  x .  f  1  ,  .  x .  , 

a  (— )  «„  +  —  b  (— )  uT 
t  e  e 


x  G  O,  u  e(x ) 


|  ffV. -)(«/)*  =  o 

2  £ 


0. 


(2.24) 


The  analysis  of  the  control  problem  involves  homogenization  of  this  system. 
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Let 


At  =  a  (y)  dyy  +  b(y)dy 


(2.25) 


with  its  formal  adjoint  defined  by 


A  i  =  dv  [  a  (y  )  dy  .  ]  -  dy  [(6  {y  )  -  ay(y  )) 


The  problem 


A  i  m  =  0,  y  — ►  m(y)  periodic 


(2.27) 


>0,  J  m(y)  dy  =  1 


has  a  unique  solution  m  (.)  on  Y  =  S°,  the  unit  circle,  with 


0  <  m  <  m(y)  <  m  <  oo. 


So  m  (.)  is  a  density  on  Y\  We  assume  that  6  (.)  is  centered 


f  m  (y  )  b  (y  )  dy  =  0. 

y 


(2.29) 


As  a  consequence  the  system 


^ix(y)=  b(y) 


(2.30) 


V  — -  X(y  )  periodic,  f  x(y)  dy  =0 
Y 


has  a  well  defined  solution,  x(')  is  the  corrector  associated  with  the  problem. 


As  before  we  look  for  u f  in  the  form 


u \x  )  —  u  l(x  ,y  )  =  u0(x  ,y  )  +  e  u  ,(x  ,y  )  +  ,  y  =  x/e 


and  we  use 


dx  <j>(x  ,y  )  =  <j>x  (ar  ,y  )  +  —  <j>y(x,y)  ,  y  =  x  /e 


dxl  <j>(x  ,y  )  =  <bJX  (x  ,y  )  +  —<f>ly{x.y)  +  <j>yy(x  ,y)  . 


Substituting  in  (2.24),  we  have 


2  1  1 

a(y)[u0„  +  —  u0xy  +  —  u0yy  ]  +  a(y)  [e  +  2  ux  +  —  u  lyy  ] 

t  t  e 

+  a(y)[e2  V  2XX  ~h  26  U  2xy  ~t~  V  2 yy  ] 

+  -j6(y)[«o*  +  j  «oy  ]  +  -i-6(y)[e«i*  +  Mij/]  (2.33) 

•f  &  (jO  [  £2  «  21  +  t  U2y  ]  -  e  [  «  o  +  £«1  +  e2  «  2  ] 

-  1-  </2(*  ,2/  )  [  («0*  +  €  «i*  +  e2  )  +  —  («0y  +  €  «iy  +  <?  « 2y  )  ]2  =  O  (eZ) 

The  last  term  is 

1  12 

-  —  92{x  ,y  )  (  (  —  i£02„  -f  —  ii0l  u0y  -f  2u0l  Mly  +  u  0Z  )  (2.34) 

c  t 

+  e(2«0l  «„  +  2«  jj  Ul!(  -f  2u  1J(  u  2y  +  2U0J  u2y)]  +  O  (e2) 

Equating  coeiTicients  of  like  powers  of  e,  we  obtain 

(«~2)  a(y)u  0yy  +  b  (y)  u0y  -  g\x  ,y)  u £y  =  0  (2.35a) 

di 

(€_1)  « (3/ )  M  lyy  +  b(y)uly  +  2a(y)u0xx  (2.35b) 

+  b  (y  )  u0x  -  g2(x  ,y)u0x  u0y  =0 

(e°)  a  (2/)  M  2j/y  ■/"  b  (y  )  ii  2y  +  “Za  (y  )  u  lxy  +  b  (y  )  u  lx  (2.35c) 

+  a(y)« 0«  -  c«o  -  ~  2(ar ,y )  u ox  -  g2{x,y)u0x  «„  =  o. 

Choosing  u0{x,y)  =  u0(x),  which  must  be  justified,  satisfies  (2.35a).  We  can  then  solve  (2.35b) 
by  choosing 

u \ix  ,y)  —  -x(y)v-oAx)  +  «,(i).  (2.36) 

Equation  (2.35c)  has  a  solution  for  u2(x  ,y  )  if 

f  ™{y)  {  -2 a  (y)xy  ^  0X2  -  b(y)xy  u  ox  +  a  (y)  ^  o xx  (2  37) 

Y 

-  cu0  -  ^  g~(x  ,y )  (  l  -  2xy  )  w02  }  dy  =0. 

This  gives  an  equation  for  u0{x  ) 
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Q  w  ojx 


c  u o  -  —  r  «0l  —  o 


(2.38) 


where 


q  =  f  m{y)  {a(y)  [l  -  2 xy{y)]  -  x(y  )  Hy )  } 
y 

r  =  /  m(y)  ?2(a:  )  [  X  -  2Xj,  (t/ )  ]  dy. 


Remark.  From  the  definition  of  A 1  and  the  corrector  x(V )  we  have 


j  m  (y  )  b  (y  )  x(V  )  dy  =  j  [  a  {y  )  xy!/  +  t>  (y)  Xy  }  X(y  )m(_y)  dy 

Y  Y 


(2.40) 


f  x(y  )  dyv  [  a  (y)  x(y)  m(y)}  dy  -  /  x(y  )  dy  [  b(y)x(y)  m(y)]  dy 

Y  Y 


Also,  using  (2.30), 


/  m  (y)  b(y)  x(y)  dy  =  f  x(y  )  a  (y  )  m  ( y  )  Xyy  (y  )  dy  (2  41) 

y  Y  ' 

-  f  x(y  )  b(y  )  m  (y  )  Xy  dy  +  2  /  x(?/ )  xy  [  «  (y  )  m  (j/ )  ]  dy 

Y  Y 

Adding  these  two  expressions,  we  have 

2  f  m(y)  b(y)x(y )  dy 

Y 

=  2  f  xiy  )  a  {y  )  m  {y  )  Xyy  dy  +2  f  X(y  )  xy  [  «  ™  ]y  dy  (2  42) 

y  y 

=  -2  /  <5j,  [  X  0  m  ]  Xj,  rfy  +  2  /  x(y )  Xy  [«»»]»  -2  f  Xy  a  (y)  m  (y  )  Xy  dy 

Y  Y  Y 

Thus,  5  may  be  rewritten  as 

<7  —  /  m  (y  )  {a  (2/ )  [  1  -  2  xy  +  Xy  ]  }  dy  (2  43) 

y 

=  /  m  (y)  a  (y)  [  l  -  Xy  i2  dy 

Y 

and  clearly  q  >  0. 

The  term  q  in  (2.39)  summarizes  the  effects  of  the  averaging  process  on  the  uncon t, rolled  sys¬ 


tem.  The  homogenization  process  interacts  with  the  control  system  through  the  term  F,  whose 
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form  would  be  difficult  to  “guess”  from  simple  averaging  procedures. 


3.  Continuum  Model  for  a  Simple  Structural  Mechanical  System 

3.1.  Problem  definition 

Consider  the  truss  shown  in  Figure  1  (undergoing  an  exaggerated  deformation).  We  shall 
assume  that  the  truss  has  a  regular  (e.g.,  triangular)  cross-section  and  no  “interlacing”  supports. 
We  assume  that  the  displacements  of  the  system  are  “small”  in  the  sense  that  no  components  in 
the  system  buckle.  We  are  interested  in  describing  the  dynamical  behavior  of  the  system  when 
the  number  of  cells  (a  unit  between  two  (triangular)  cross  sections)  is  large;  that  is,  in  the  limit  as 

(  =  l  /L  —*  0.  (3.1) 


Figure  1.  Deformed  truss  with  regular  cross-section. 
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We  shall  make  several  assumptions  to  simplify  the  analysis.  First,  we  shall  assume  that  the 
triangular  sections  are  essentially  rigid,  and  that  all  mobility  of  the  system  derives  from  the  flexi¬ 
bility  of  the  members  connecting  the  triangular  components.  Second,  we  shall  ignore  damping 
and  frictional  effects  in  the  system.  Third,  we  shall  confine  attention  to  small  transverse  displace¬ 
ments  and  small  in  plane  rotations  <j>{t  ,x)  as  indicated  in  Figure  1,  ignoring  longitudinal 

and  out  of  plane  motions  and  torsional  twisting.  Fourth,  we  shall  assume  that  the  mass  of  the 
triangular  cross  members  dominates  the  mass  of  the  interconnecting  links. 

Systems  of  this  type  have  been  considered  in  several  papers  including  (Noor  et  al.  1978) 
(Nayfeh  and  Hefzy  1978)  (Anderson  1981)  and  (Renton  1984).  In  those  papers  a  continuum  beam 
model  was  hypothesized  and  effective  values  for  the  continuum  system  parameters  were  computed 
by  averaging  the  associated  parameters  of  the  discrete  system.  Our  approach  to  the  problem  is 
based  on  homogenization-asymptotic  analysis  and  is  quite  different. 

The  assumptions  simplify  the  problem  substantially,  by  suppressing  the  geometric  structure 
of  the  truss.  We  can  retain  this  structure  by  writing  dynamical  equations  for  the  nodal  displace¬ 
ments  of  the  truss  members.  For  triangular  cross  sections  nine  parameters  describe  the  displace¬ 
ments  of  each  sectional  element.  The  analysis  which  follows  may  be  carried  over  to  this  case,  but 
the  algebraic  complexity  prevents  a  clear  presentation  of  the  main  ideas.  As  suggested  in  (Noor 
et  al.  1978)  one  should  use  a  symbolic  manipulation  program  like  MACSYMA  to  carry  out  the 
complete  details  of  the  calculations.  We  shall  take  up  this  problem  on  another  occasion;  for  now 
we  shall  treat  the  highly  simplified  problem  which,  as  we  shall  see,  leads  to  the  one  dimensional 
Timoshenko  beam. 

We  shall  begin  by  reformulating  the  system  in  terms  of  a  discrete  element  model  as  sug¬ 
gested  in  (Crandal  et  al.  1980);  see  Figure  2.  In  this  model  we  follow  the  displacement  ;?,  (/  )  and 
rotation  <^,  (2)  of  the  ith  mass  M.  The  bending  springs  (/q*)  tend  to  keep  the  system  straight  by 
keeping  the  masses  parallel  and  the  shearing  springs  (k'e)  tend  to  keep  the  masses  perpendicular  to 
the  connecting  links.  We  assume  small  displacements  and  rotations  so  the  approximations 

sin  )  (3.2) 
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Figure  2.  A  lumped  parameter  model  of  the  simplified  truss  system 


tan  1  [>?,•(<  )A  ]  =  Vi(t)/I 


are  valid. 


In  this  case  the  (approximate)  equations  of  motion  of  the  ith  mass  are0 


d~4>i  ~  1  ,i(l  +  )  -  rh  (t)  , 

— —  =  —  K  {  1  - 7 - J  -  )  \ 

dt~  r  l 

+  -  (j);  (t) 

+  sr  { ki  [  *'+iK  >  {  !.  ] } 

f/2?h  f  „■  ,  »7,(0  -  ,  j,,,,, 

— 7- —  =  [ - 7 - ]  ^  &  (t )  ) 

at  -  1 

where  we  have  normalized  M  =  1  and  defined 


(3.3b) 


°The  spring  constants  depend  on  l  since  they  represent  the  restorative  forces  of  flexed  bars,  bent  by  different 
amounts.) 
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(3.4) 


m  ± 


[  Vi-i  -  Vi 


and  similarly  for  Si  4>i  ■ 

To  proceed,  we  shall  introduce  the  nondimensional  variable  t  —  l  /L  and  rewrite  the  sys¬ 
tem  (3.3)  as 


r 


dt2 


=  y  Ki  {  V£  +  ViV  )  -  4>i  }  +  v£  +  { Ki  vf  )  } 
=  -  V"  {  Ki  [  v(  +  v/U )  -  Mt )  ] } 


(3.5) 


where 

Ki  =  /c/  L  ,  A7  =  L  (3.6) 

v£  +  Vi  =  —  [  Vi+ 1  -  >?,•  ]  .  Vc“  Vi  =  —  [  Vi  -  Vi-1  ]  • 
e  e 

Normalizing  L  —  1,  we  associate  a  position  2  €E  (  -  — ,  —  |  with  each  mass;  and  we  introduce 

2  2 

the  notation 


i){t  ,x{)  =  Vi(t),  4>(t  ,x{)  =  4>i{t )  .  (3.7) 

Having  normalized  L  =  1,  we  have  e  =  /  and  x,-+1  =  2,-  +  l  —  2,-  +  e.  Let  z?  =  {2,}  be 

the  set  of  all  points  in  the  system.  In  this  notation 


(ve  +  v)((  >2  )  =  —  [  v(t  ,2  +  <0  -  v(t  .2  )  ] 

(Ve~  »7)(<.2  )  =  —  [  »;(<  ,2)  -  »/(<  ,x  -  t)  ],  x  €  Z 
e 


(3.8) 


and  the  system  is 


rf20f(f  ,2,-  ) 


r//2 


d2Vl{t  ,2,-) 


/C  (2,-  )  {  Vf  +  *?'(<  .2,-  )  -  4>C(t.  ,X{  )  } 

+  r  vt+  {  /v*  (2.  )  Sfc~  4>({t.  ,x{)  } 
vf  {/7  (2,  )  [  v{  +  >yf(<  >2,- )  -  4>€{t  ,2,-)]},  2  e  £ 


(3.8) 
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The  scaling  of  (3.9)  may  be  interpreted  in  the  following  way:  Formally,  at  least,  the  right 
sides  of  both  terms  in  (3.9)  are  0  (e-2).  This  implies  that  the  time  variations  are  taking  place  in 
the  “fast  time  scale”  r  —  t  /e.  Also,  the  spatial  variations  are  taking  place  in  the  “microscopic 
scale”  x  which  varies  in  e-increments  (e.g.,  z,-+i  =  a:,  +  e).  Introducing  the  macroscopic  scale 

z  —  e  x  ,  and  the  slow  time  scale  a  —  e  t,  we  may  rescale  (3.9)  and  observe  its  dynamical  evolu¬ 
tion  on  the  large  space-time  scale  on  which  macroscopic  events  (e.g.,  “distributed  phenomena”) 
take  place. 

Rewritten  in  this  spatial  scale,  the  system  becomes 


e 

1? 


±K,(—){6(+  r^t,—)  -  <t>V,~ )} 


+  ±r+{r  Kb(—)8(-tV,—)} 


(3.10a) 


_ € 

dt 2 


=  4  8‘- 

c~ 


{K,(-±)[8t  +  -  etV,-^-)}} 


(3.10b) 


where 


6ci  =  €  ve±  —  O  (1)  in  e. 


(3.11) 


The  essential  mathematical  problem  is  to  analyze  the  solutions  <j)e ,  »7£  of  (3.10)  in  the  limit  as 
e  — ►  0. 


3.2.  Mathematical  analysis 

To  proceed,  we  shall  generalize  the  problem  (3.10)  slightly  by  allowing  Ke  and  Kb  to 
depend  on  z  as  well  as  z  /e.  This  permits  the  restoring  forces  in  the  model  system  to  depend  on 
the  large  scale  shape  of  the  structure  as  well  as  on  local  deformations.  We  use  the  method  of 
multiple  scales;  that  is,  we  look  for  solutions  of  (3.10)  in  the  form 

rj({t)  =  if{t  ,z  ,y),  </>e(t)  =  <j>V  ,z  ,y),  y  =  (3.12) 
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and  we  have 


K,  —  K,  (z  ,y ),  Kb  =Kb(z,y),  y  =  j-  (3.13) 

On  smooth  functions  ib(z  ,—)  the  operators  8C±  satisfy 

e 

(8C  +  ip)(z  ,y  )  =  if(z  +e,y  +1)  -  ip(z  ,y  ) 

—  0(2  ,y  +1)  -  ip{z  ,y)  +  if{z  +  e,  y  + 1)  -  if(z,y+l)  (3.14a) 

==  (-S'  +  ip)(z,y)  +  e~{z,y+\)  +  -i-  e"^-^{z  ,y  +1)  +  0  (e3) 
az  2  dz2 

(t>c~  4>)(z  ,y)  =  4>{z  ,y)  -  ip(z-e,y- 1) 

=  ip(z  ,y )  -  ip(z,y- 1)  +  ip{z  ,  y-l)  -  tp{z-e,y- 1) 

=  (51 "  f>){z  ,y  )  -  e^-{z,y+l)  +  ^-t2^~{z,y- 1)  +  0  (e3)  (3.14b) 

az  2  dz 

We  assume  that  0£  and  if  may  be  represented  by 

,z  ,y)  =  <f>0(t  ,z)  +  e  f>i(t  ,z  ,y  )  +  ■■■  (3.15) 

if(t  ,z  ,y)  =  r)0{t  ,z)  +  t  »;,(/  ,z  ,y)  + 

and  substituting  (3.15)  in  (3.10)  and  using  (3.13)  (3.14),  we  arrive  at  a  sequence  of  equations  for 
(0O.  Vo)>  (0i,  Vi),  •  -  -  by  equating  tlie  coefficients  of  like  powers  of  e. 

Starting  with  e-2,  e_1,  e°,  •  •  •  ,  we  have 

S  +  [r  K{z  ,y)  S~  0O(<  ,z )  ]  =  0  (3.16) 

which  is  trivially  true  from  (3.14b)  (3.15).  The  same  term  involving  t/0(f  ,x  )  from  (3.10b)  is  trivi¬ 
ally  satisfied  by  the  assumption  (3.15).  Continuing 

■j  l  S  +  {  rICb  (z  ,y)S~  4>ft  ,z  ,y)}  (3.17) 

+  Ke(z  ,y)  {  S  +  r)0(t  ,z  )  -  <j>0(t  ,z  )  }  )  =  O 
which  may  be  solved  by  using  the  corrector  x<f>(z  ’V  )  and  taking 
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Mt.z.y)  =  xAz  ’V )  0o(<  >z) 


(3.18) 


with 

S  +  {rKb(z  ,y)  S  ~x^.z  ,y)  }  =  K,(z  ,y)  (3.19) 

If  we  regard  2  as  a  parameter  in  (3.19),  then  there  exists  a  solution  x</»  unique  up  to  an  additive 
constant,  if  Kb  (z  ,  .  ),  K,  (z , .  )  are  periodic  in  y,  if  there  exist  constants  A  and  B  so  that 

0  <  A  <  Kb(z  ,y)  <  B  <  oo  (3.20) 

and  if  the  average  of  K„  (z  ,  .  )  is  zero 

L/Z 

—  I  K.(s  ,y)  dy  =  0  (3.21) 

h  -L/2 

which  holds  if  the  system  is  pinned  at  the  ends  as  indicated  in  Figure  2.  Let  us  assume  that 
(3.20)  (3.21)  hold,  and 

0  <  A  <  I<,(z  ,y)  <  B  <  oo  (3.22) 

(which  we  shall  need  shortly). 

Considering  (3.10b),  the  O  (e"1)  term  in  the  asymptotic  expansion  is 

~[S-  {K,{z,y){S+  Ih{t,z,y)  -  &,(<,*))}]=  0  .  (3.23) 

Again  we  introduce  the  corrector  x<iiz  ’V  )>  anc*  take  ^  in  the  form 

hit  >z  ,y )  =  x„(z  >y )  ,z )  (3.24) 

which  gives  the  equation  for  the  corrector 

S  ~  {  Ke{z  ,y)  \  S  +  x„(z  ,V)  -  1  ]  }  =  0  (3.25) 

or 

S  -  {  K„{z,y)  S  +  x,,(z,y)  }  =  Ke{z,y)  -  Ke(z,y- 1)  (3.20) 

By  hypothesis  the  right  side  in  (3.26)  is  periodic  in  y  and  has  zero  average  (3.21).  Hence,  (3.20) 
has  a  periodic  solution,  unique  to  an  additive  constant. 
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Continuing,  the  O  (e°)  term  in  (3.10a)  is 


5  +  {rKb  (z  ,y)  S  ~  <f>2{t  ,z  ,y )  }  +  K,  (z  ,y  )  [  S  +  f?,(<  ,z,y)  -  ^(t  ,z,y)] 

+  K,  (z  ,y  )  —£—{t  ,z  )  +  S+  {rKb  (z  ,y  )  -J— ^(f  ,z  ,y  )  }  (3.27) 

OZ  oz 

+  S  +  {rKb(z  ,y)  ~^4>0{t  ,z  ,y)}  +  —  {rKb  (z  ,y  +1)  }  —  <p0{t  ,z  ) 

d2  d2<j>0 

+  -r-r  {  rKb  {z  ,y  +1)  }  M*  >z  )  -  =  0  ■ 

dz 2  at* 

This  should  be  regarded  as  an  equation  for  <j>2  as  a  function  of  y  with  (l,z)  as  parameters.  In  this 
sense  the  solvability  condition  is  as  before,  the  average  of  the  sum  of  all  terms  on  the  left  in 
(3.27),  except  the  first,  should  be  zero.  We  must  choose  <}>0  so  that  this  in  fact  occurs;  and  that 
defines  the  limiting  system. 

Using  the  correctors  (3.18)  (3.24),  we  must  have 


d260  d2(j)0 

Average  {  -  -  [  S  +  (rKb  (z  ,y ))  +  S  (rKb  (z  ,y )  \^(z  ,y  ))  ] 

(</)  dt2  dz~ 


d<b0  d  dt]0 

~d7  ^  ~dz  <"rKb  ^  ,V  +1^  ^  “  ~d7K‘<<Z’y'> 


(3.28) 


-  4>0  l  ~~H  (rKt>  (z  >V  +1))  +  s  +  (rl<b(z,y)  -^x4z’V)) 


dz‘ 


+  K,  (z  ,y  )  (S  +  x,t(z  ’V  )  -  x4z  -y  ))  ]  }  =  0 


Defining  the  functions  El  ( z  ), 

d2<j>0  _ 
dt 2 


G  (z  )  by  the  associated  averages  in  (3.28),  the  averaged  equation  is 

JL(EI(z)?p-)  +  G(z )—  -  o  (3.29) 

dz  oz  oz 


which  is  the  angular  component  of  the  Timoshenko  beam  system  (Crandall  et.  al.  1980  p.  348). 

Arguing  in  a  similar  fashion,  we  can  derive  the  equation  for  the  macroscopic  approximation 
displacement  of  the  lattice  system  in  terms  of  the  “equivalent”  displacement  r/0{t,z)  in  the 
Timoshenko  beam  system 


d2>?o  __  _<9_ 
dt2  dz 


[  N  {z  )  ( 


dri0 

dz 


Mt  >z)  )  1 


(3.30) 
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3.3.  Summary 


We  have  shown  that  a  simplified  model  of  the  dynamics  of  the  truss  with  rigid  cross  sec¬ 
tional  area  may  be  well  approximated  by  the  Timoshenko  beam  model  in  the  limit  as  the  number 
of  cells  (proportional  to  L  /l)  becomes  large.  The  continuum  beam  model  emerges  naturally  in 
the  analysis,  as  a  consequence  of  the  periodicity  and  the  scaling. 

To  compute  the  approximate  continuum  model,  one  must  solve  (3.19)  and  (3.26)  (numeri¬ 
cally)  for  the  correctors  and  then  compute  the  parameters  in  (3.29)  (3.30)  by  numerically  averag¬ 
ing  the  quantities  in  (3.28)  (and  its  analog  for  (3.10b))  which  involve  the  correctors  and  the  data 
of  the  problem. 

4.  Homogenization  and  Stabilizing  Control  of  Lattice  Structures 

In  this  section  we- show  that  the  process  of  deriving  effective  “continuum”  approximations 
to  complex  systems  may  be  developed  in  the  context  of  optimal  control  designs  for  those  systems. 
This  procedure  is  more  effective  than  the  procedure  of  first  deriving  homogeneous/continuum 
approximations  for  the  structure,  designing  a  control  algorithm  for  the  idealized  model,  and  then 
adapting  the  algorithm  to  the  physical  model.  In  fact,  separation  of  optimization  and  asymptotic 
analysis  can  lead  to  incorrect  algorithms  or  ineffective  approximations,  particularly  in  control 
problems  where  nonlinear  analysis  (e.g.,  of  the  Bellman  dynamic  programming  equation)  is 
required. 

We  shall  apply  the  combined  homogenization  -  optimization  procedure  described  in  section 
1  (based  on  (Bensoussan,  Boccardo,  and  Murat  1984))  to  the  problem  of  controlling  the  dynamics 
of  lattice  structures  like  the  truss  structure  analyzed  in  the  previous  section.  We  shall  only  for¬ 
mulate  a  prototype  problem  of  this  type  and  discuss  its  essential  features. 

Consider  the  model  for  the  lattice  structure  analyzed  in  section  3  with  control  actuators 
added.  The  truss  shown  in  Figure  1  is  again  constrained  to  move  in  the  plane  and  torsional 
motion  is  excluded  to  simplify  the  model  and  confine  attention  to  the  basic  ideas.  Now,  however, 
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we  include  a  finite  number  of  actuators  acting  to  cause  transverse  motions.  The  truss  with  actua¬ 
tor  forces  indicated  by  arrows  is  shown  in  Figure  3.  The  corresponding  discrete  element  model  is 
shown  in  Figure  4. 


Suppose  that  the  physical  actuators  act  along  the  local  normal  to  the  truss  midline  as  shown 
in  the  figures,  and  that  the  forces  are  small  so  that  linear  approximations  to  transcendental  func¬ 
tions  (e.g.,  sin  <j>{  =  <f>i ,  etc.)  are  valid.  Then  the  controlled  equations  of  motion  of  the  discrete 
element  system  are  (recall  equation  (3.5)) 


dt 


<t>;  l 

=  -k:  {ve 


ruV)  -  *n  +  vf+  {Ki  d>iV)  } 


(4.1) 


Figure  3.  Truss  with  transverse  actuator  forces. 
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Figure  4.  Discrete  element  model  of  tlie  controlled  truss. 


dt 2 


-  Ve "  {  Ki  [  V£  +  »?/«  )  -  4>Rt  )]}+"£  )  *s(t ) 

j  =1 


where  the  notation  in  (3.6)  has  been  used. 


J ) 


O  iVi 

x  i  =  y 


(4.2) 


and  *  j  ,  j  —  are  the  locations  of  the  actuators.  Ilence,  if-£(t",»’y)  =  0  for  all  /  =  1 

there  is  no  actuator  located  at  the  ith  point  which  corresponds  to  the  physical  point  x  E  [0 ,L}. 
The  number  m  of  actuators  is  given  at  the  outset  and  does  not,  of  course,  vary  with  the  scaling. 

The  control  problem  is  to  select  the  actuator  forces  as  functions  of  the  displacements  and 
velocities  of  components  of  the  structure  to  damp  out  motions  of  the  structure.  Measurements 
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would  typically  be  available  from  a  finite  number  of  sensors  located  along  the  structure.  We  shall 
not  elaborate  on  this  component  of  the  model,  and  shall  instead  assume  that  the  entire  state  can 
be  measured.  To  achieve  the  stabilization,  we  shall  associate  a  cost  functional  with  the  system 
(4.1).  Let 

«(0  =  I  «l(0.  •••>  «m(f)]r  (4.3) 

be  the  vector  of  control  forces,  and 

00  N 

AM  «(•)]  =  /  £  { «.  1 )  ?  +  Mtf(0]2 

0  1=1 

+  <*.•  [  4><(t )  ]2  +  A  [  mV )  ]2  (4.4) 

m 

+  Yj  ^(*  ,ij  )  M 2(f  )  }  e  dt 

)  =i 

where  (a,-,  6.  )  and  (a,-,  A)  are  non-negative  weights.  Formally,  the  control  problem  is  to  select 
8(i  ,ij)  Uj(t),  i  —  ,  j  —  1  ,...,m  to  achieve 

Tj)  r  [  U  (,)  5  (4.5) 

subject  to  (4.1)  (4.2)  and  the  appropriate  boundary  conditions.  The  case  7  -+  0  corresponds  to 
stabilization  by  feedback. 

The  analysis  of  this  control  problem  is  based  on  the  scaling  used  in  section  3,  equations  (3.5) 

-  (3.11).  Let  r  —  t  /t  be  the  fast  time  scale,  then 

00  N 

J1  [«(•)]  =  /  e  £  {  «•-  [  M'M  ]2  +  b;  [  fj/(r)  ]2 
0  »  =1 

+  a,-  e2lhc(r)f  +  A  e2[  ]2  (4.6) 

m 

+  8{i,ij)v-{T)}  (I  ~ nT  d  T. 

y=i 

with  0/(r)  ==  <j>;(cr),  etc. 

Let  {(j),  j)  be  the  state  vector  of  the  system  (4.1)  with  <j>  =  [  <j>lt  •  •  ■  ,  <j>N  }T  and  simi¬ 
larly  for  the  other  terms.  Let  V  =  Vtn((j),  be  the  optimal  value  function  for  the  problem 

(4.1)  (4.6).  Then  the  Bellman  problem  associated  with  (4.1)  (4.6)  is 
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(4.7) 


N 


«  E  \<t>vK  +  m  v%  } 


N 


+  €  s  {  K'  ( v£  +  m  -  A  )  +  -  v£  M  *V  v£“  4>i  ] }  vf 


N 


+  f  E  {  -v£_  [^»‘(v£  *?,•  -  A)1  }E), 


iV  m 


mi^  {eE  E  t^(*  .*>)«;•  +  <*(»  .*;•)“/]  } 


*j  G  ,-=l  y=i 


JV 


+  £  E  [  «.■  A2  +  h  2  +  f2  ( «,•  ^2  +  A  »?,-2 )  ]  -  =o. 


REMARKS: 

(1)  Note  that  the  minimization  in  (4.7)  is  well  defined  if  the  admissible  range  of  the  control 
forces  is  convex  since  the  performance  measure  has  been  assumed  to  be  quadratic  in  the  con¬ 
trol  variables  6(i  ,:'y  )  My  . 

(2)  Since  we  have  not  included  the  effects  of  noise  in  the  model,  the  state  equations  are  deter¬ 
ministic  and  the  Bellman  equation  (4.7)  is  a  first  order  system.  To  “regularize”  the 
analysis,  at  least  along  the  lines  followed  in  conventional  homogenization  analysis,  it  is  use¬ 
ful  to  include  the  effects  of  noise  in  the  model  and  exploit  the  resulting  coercivit3r  properties 
in  the  asymptotic  analysis. 

(3)  If  we  introduce  the  macroscopic  spatial  scale  z  =  ex  ,  the  mesh  {rr,-  },  and  the  variables 

<t>{t  ,Zj )  =  ),  4>{t  ,z{ )  =  4>i\t ),  etc.  (4.8) 

then  the  sums  may  be  regarded  as  Riemann  approximations  to  integrals  over  the  macros¬ 
copic  spatial  scale  The  asymptotic  analysis  of  (4.7)  with  this  interpretation  defines  the 
mathematical  problem  constituting  simultaneous  homogenization  -  optimization  for  this 
case. 
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5.  Effective  conductivity  of  a  periodic  lattice 


In  this  section  we  consider  a  version  of  a  heat  conduction  problem  treated  by  Kunnemann. 
Simple  expressions  for  thermal  properties  of  composite  materials,  have  been  derived  in  the  past 
using  homogenization  techniques.  The  derivation  of  effective  conductivities  for  discrete  structures 
is  useful  for  assessing  the  behavior  of  such  structures  in  variable  environmental  conditions. 

5.1.  Problem  definition 

Let  Z  —  {o,  ±l,  ±2,  •  •  •  }  and  Zd  —  Z  X  •  •  •  X  Z  ( d  times)  be  a  (/-dimensional  lat- 
tice.  Let  e  >  0  be  a  number. small  relative  to  1.  We  want  to  describe  the  effective  conduction  of 
thermal  energy  on  the  e-spaced  lattice  tZd  .  Let  t{  =  (0,0, ...0,1,0, ...0)r  with  1  in  the  ith  posi¬ 
tion,  i  —  1,2, ...,</.  If  x  is  a  point  in  eZd ,  then  x  ±  ee,-,  1  <  i  <  d ,  are  the  nearest  neighbors  of 
x.  Let  a  ±(x ),  x  EZd,l<i<d,be  the  two  functions  defined  on  the  lattice,  and  assume 

«,•  (x  )  =  a  +(x)  —  a,-_(x  +  e,  ),  x  6  Zd ,  1  <  j  <  d  (5.1a) 

0  <  A  <  a{(x )  <  B  <  oo,  \-/  x  £  Zd ,  1  <  i  <  d  (5.1b) 

a,  (x)  is  periodic  with  period  l  >  1  (5.1c) 

in  each  direction,7  1  <  i  <  d. 

Next  let 

ai±(x )—  ai±(~),  x  &  tZd ,  1  <  i  <  d.  (5.2) 

Equation  (5.1b)  means  that  the  conduction  process  is  reversible  and  that  the  conductivity 
a,  (x)  is  a  “bond  conductivity,”  i.e.,  independent  of  the  direction  in  which  the  bond  (x  ,x  +  e,- ) 
is  used  by  the  process.  Equation  (5.2)  means  that  the  configuration  of  bond  conductivities  a{e±  (.) 
on  e  Zd  is  simply  a,-±(.)  on  e  Zd  “viewed  from  a  distance.”  Assumption  (5.1c)  imposes  a  regular¬ 
ity  condition  on  the  physics  of  the  conduction  process.  An  assumption  like  this  is  essential  for 
existence  of  a  limit  as  e  — +  0.  In  one  dimension  the  situation  is  illustrated  in  Figure  5a, b.  A  sys¬ 
tem  similar  to  this  with  random  bond  conductivities  was  treated  by  Kunnemann  (1983)  by 
7The  period  may  be  different  in  different  directions. 
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imposing  some  ergodicity  properties  on  the  bond  conductivities. 

One  can  associate  with  this  system  a  random  (jump)  process 

{  Xc(t  ,x ),  t  >  0,  x  G  eZd  } 

on  the  €-spaced  lattice.8  In  effect,  as  e  — *•  0  ,  {  Xe  }  converges  to  a  Brownian  motion  on  the  lat- 
tice;  and  the  main  result  of  the  analysis  is  an  expression  for  the  diffusion  matrix 
Q  =  [  9,y ;  t  ,j  =  1,2 ]  of  this  process.  This  matrix  describes  the  macroscopic  diffusion  of 


^eflniUon  of  this  process  Is  not  necessary  for  the  analysis,  but  it  bolsters  the  Intuition. 
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thermal  energy  in  the  system.  It  is  the  effective  conductivity. 

We  shall  carry  out  the  asymptotic  analysis  of  this  system  in  the  limit  as  e  — *  0  using  homo¬ 
genization.  Let 


(V/ 
(V/  + 


for  any  u  square  summable  on  e 
tor  in  IR^  .  Then 


u  )(x )  =  —  [  u  (x  -  ee,  )  - 
e 

u  )(x)  =  —  [  « (x  +  ee,-)  - 
e 

a;  £  e  Zd  ,  1  <  «  <  d , 

Zd  or  square  integrable  on  IR^ 


«  (*)  ] 

«  (z)  ] 

with  e,-  the  ith 


(5.3a) 


(5.3b) 


natural  basis  vec- 


du  c(t  ,x  ) 

at 


S  V/“  [  «.-(—)  V,£+  u\t,x)} 
i  =  i  e 

t  L£  «£(f  ,z) 


(5.4) 


is  the  diffusion  equation  on  the  e-spaced  lattice  with  density  u\t  ,x)  and  conductivity  a,-(x/(). 
We  are  interested  in  an  effective  parameter  representation  of  the  thermal  conduction  process  as 
c  -+  0. 

Remark:  Although  probabilistic  methods  are  not  required  in  the  analysis,  the  associated  proba¬ 
bilistic  framework  has  a  great  deal  of  intuitive  appeal.  The  operator  L(  may  be  identified  as  the 
infinitesimal  generator  of  a  pure  jump  process  Xe(s  )  in  the  “slow”  time  scale  s  =  e2  t ;  (Brei- 
man  1968).  Moreover,  L£  is  selfadjoint  on  e  Zd  with  the  inner  product 

(/.*)  =  £  /(*)*(*)•  (5.5) 

x  £Zd 

Hence,  the  backward  and  forward  equations  for  the  process  X€(s  )  are,  respectively, 

dPc(y>U.xl  =  [L£P£(y,f  [  .)](*)  (5.6) 

at 

dp((y’l  I  =  [L‘pe(.,t  |*)](y) 

at 

So  the  process  is  “symmetric”  in  the  sense  of  Markov  processes  (Breiman  1968). 
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The  asymptotic  analysis  of  (5.4),  when  interpreted  in  this  context,  means  that  as  the  bond 
lattice  is  contracted  by  c  and  time  is  sped  up  by  e~2,  the  jump  process  {  Xe(s  )  }  approaches  a 
diffusion  process  with  diffusion  matrix  Q  .  In  other  words,  on  the  microscopic  scale  thermal 
energy  is  transmitted  through  the  lattice  by  a  jump  process;  but  when  viewed  on  a  macroscopic 
scale  the  energy  appears  to  diffuse  throughout  the  lattice.  The  microscopic  physics  are  described 
in  (Kirkpatrick  1973)  and  (ICittel  1976).  The  approximation  developed  below  for  a  periodic  lattice 
is  similar  to  the  one  developed  by  Kunnemann  for  a  random  lattice.  This  similarity  demonstrates 
the  robustness  of  the  method,  and  the  limted  dependence  of  the  macroscopic  properties  of  the 
medium  on  the  details  of  the  microscopic  variations  of  the  structure. 

Because  the  basic  problem  (5.4)  is  “parabolic,”  we  can  introduce  the  probabilistic  mechan¬ 
ism  and  make  use  of  it  in  the  analysis.  In  the  “hyperbolic,”  structural  mechanical  problems  we 
treated  before  this  device  is  not  available. 


5.2.  Asymptotic  analysis-homogenization 


The  essential  mathematical  step  is  to  show  strong  convergence  of  the  semigroup  of  Le,  say 


A  „Lt 


T\t)  =  e - ►  T(t)  t  e 

e->0 


(5.7) 


and  to  identify  the  limiting  operator 


*  ..  d2 

,  j  =  i  9'3  dx<  dxi 


(5.8) 


This  is  accomplished  by  proving  convergence  of  the  resolvents 


for  a  >  0,  [  -  Lc  +  a 


-L  + 


(5.9) 


That  is,  if  /  is  a  given  function  and 


n€{.)  t  [  -Lc  +  a}-1  f 
«  (  •  )  =  [  ~L  +  a)-1  f 


(5.10) 


then  v  —*  u  in  an  appropiate  sense. 
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The  method  of  multiple  scales  will  be  used  to  compute  the  limit.  Because  the  conductivities 
a,  (z)  in  (5.4)  do  not  depend  on  time,  we  may  work  directly  with  L(  rather  than  the  parabolic 
PDE  (5.4)  (cf.  (Bensoussan,  Lions,  and  Papanicolaou  1978)  Remark  1.6,  p.  242).  The  method  of 
multiple  scales  is  convenient  because  it  is  a  systematic  way  of  arriving  at  the  “right  answers”  - 
something  which  is  not  always  simple  in  this  analysis. 

Bearing  in  mind  (5.10),  we  consider 

(Lf  «*)(*)  =  /  (*)  (5.11) 

with  uc(x  )  in  the  form 

u  e(x  )  ==  u0(x  ,—)  +  e  u  i(z  ,  — )  +  €2  u  2(x  ,—)  ~b  ■  ■  ■  (5.12) 

e  t  e 

with  the  functions  Uj{x,y)  periodic  in  y  E  e  Zd  for  every  j  =  0,1 .  (As  it  turns  out  the 

boundary  conditions  are  somewhat  irrelevant  to  the  construction  of  “right  answers.”)  To  present 
the  computations  in  a  simple  form,  it  is  convenient  to  introduce  y  —  x  /(..  to  treat  x  and  y  as 
independent  variables,  and  to  replace  y  by  i/(  at  the  end. 

Recall  the  operators  v«e±  Rom  (5-4)-  Applied  to  a  smooth  function  v  —  u  (x  r/e),  we 

have 


(  v/_  u  )(x  ,y)  —  —\u(x  -  eet- ,  y  -  e,- )  -  u  {x  ,y  )  ] 


(5.13) 


=  —  [u{x  ,y  -  e,- )  -  u  (x  ,y  )  ]  +  —  [  u  (x  -  ee{ ,  y  -  e,- )  -  u  (x  ,y  -  et- )  ] 


=  —  (  V<  «  ,y  ) 
e 


du 


1  dzu 


-  c‘}  +  e  2 


{x  ,y  -  e,  )  +  0  (e2) 


where  on  functions  <f>  —  4>{y  ) 


(V,  <i>)(y)  =  <Ky  -  e.)  -  (j){y) 


(5.14) 


Defining 

(v,-+  4>){y )  =  <t>(y  +  e.)  -  4>{y)  (5.15) 

we  also  have 
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(v.e  +  u)(x,y)  =  J  (V.+  «)(*.y)  +  <y  +  e>)  (5.16) 

+  eTT^7^,2/  +  ei)  +  °(f2)- 

2  dx,-2 

Now  we  substitute  (5.12)  into  (5.11)  and  use  the  rules  (5.14)  (5.15).  Equating  coefficients  of 
like  powers  of  e,  this  leads  to  a  sequence  of  equations  for  u0,  vlr  •  •  •  .  Specifically,  (using  the 
summation  convention) 

(4£  u c)(x,y)  =  -  v<f~  !  «,-(y)  V,£+  «f] 

1  du  o 

.  ==  -7  vr  [  «,(y)  v,+  « 0(*.y) ]-  v.£~  [  «,( if )  -r — (x >y  +  e.)] 

e2  axi 

it 

-  4€V,-£“  [  a,-(y)  — |-(*.y  +  ei  )  1  +  O(e)  (5.17) 

2  3x,- 

-  —  vr  1  «.  (y  )  V,+  «i(*  ,2f )  j 

e 

du ! 

-  e  v,£  [  «,-(y )  >y  +  e>-)i  +  0 

-  vr  1  «,(</)  V,+  u2(x  ,y)}  +  0  (c)  —  f  (x) 

That  is.  labeling  each  term  by  its  order  in  e 

(r2)  Vf“  [«.’(»)  Vf+  «o)  =  0  (5.18) 

du0 

(e  ')  «  V.£  [  a.'  (y  )  - (x ’V  +  e>')l  +  V.'  [«,'(»)Vi  «i  (x,y)}=0  (5.19) 

OXi 

and  (recall  e  V.£±  is  O  (1)  in  e) 

„  i  d2u0  du  i 

(e  )  —  e  V.£  [  «,-  ( y  )  — — r(x  >V  +  «.)]  -  «  Vf  “  [  «» ( 2/ )  — (*  >V  +  e< )]  (5.20) 

2  dx2  OXi 

-  vr  [  a.(j/)  V.'+«2(:r^)]  =  /  (*) 

From  (5.18)  we  have 

a{(y  -  ei)[u0(x,y)  -  u0{x,y  -  e,- )  ]  (5.21) 

-  a;(y)[u0(x,y  +  e, )  -  u0(x,y))  =  0 

If  we  take  u0(x  ,y  )  =  u0(x  ),  this  is  trivially  true;  and  (5.19)  simplifies  to 
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(5.22) 


e  Vi 


«.  {y  ) 


3u  c 
3  a:,- 


-(*)]  +  Vi  [  a((y  )  V,+  ux{x,y) ]  =  0  . 


At  this  point  we  introduce  “correctors.”  That  is,  we  assume 

du0 

«i(z.SO=  E  Xk(y)^ -  +  "iOO  (5-23) 

*  =i  d** 

with  X*  (  •  )  the  correctors.  Using  this  in  (5.22),  we  have  (again  using  the  summation  convention) 

Vf'  [  «,  (3/ )  Vi  +  X*(y)]  4—^  +  [  a*(y  -  ek)  -  ak(y)}  =0  (5.24) 

ox k  dxk 


If  we  take  Xk(y)  35  the  solution  of 

v."  l  My)  v,+  x*(y)]  +  [«t-(y  -  «*)  -  (y ) )  ==  o  (5.25) 

(we  have  to  verify  the  well-posedness  of  (5.25)),  then  (5.24)  is  satisfied.  (The  term  u  ^rr)  is  deter¬ 
mined  (formally)  from  the  0(e)  term  in  the  system  (5.12)  (5.17).) 

Regarding  the  well-posedness  of  (5.25),  note  that 

Vi"  [  «.■  (y )  v>  +  4>(y )  ]  =  4>{y )  (5-26) 


has  a  periodic  solution  on  f  Z  which  is  unique  up  to  an  additive  constant  iff  the  average  of  the 
function  ip(y )  over  a  period  (el)  is  zero;  i.e., 


v-  £ 


l 

E 

k  =1 


V>(y  +  ken) 


1,2 


(5.27) 


This  condition  clearly  holds  in  (5.25),  and  so,  Xk(y)  *s  well  defined  (up  to  an  additive  constant). 

We  shall  determine  the  equation  for  v  0(x )  by  using  (5.23)  (5.25)  in  (5.20).  Using  the 
Kronecker  delta  function  <5,>  ,  we  have 

i  cfup  d2u, 


—  e  V/  !  «i(y  )  fa  )  .. 

2  ClX;  dxk 


-  e  f  Vi  [  ai  (y)  Xkiy  +  et- )  j 


=  Vi  [  «.  (y)  5,-*  ]  -  V,-  1  «,-(y )  Vi+Xk  \  } 


dx;  Dxi. 
d2ua 


=  /  (*) 


3a;,-  Sa:^ 


(5.28) 


^2  y 

Vi"  l«i(y)  x*(y)  ]  a  - Vi"  l  a,  (y)  Vi+  «»]  =  /  (*) 

3x,-  d.rh 


30 


The  term  in  braces  is  zero  from  (5.25).  To  obtain  the  solvability  condition  (5.27)  for  u2  in  (5.28), 


we  introduce  the  average 


£  symmetric  part  {  -  Vi  [  a.'  (y  )  Xk  (y  )  1} 


(5.29) 


Then  solvability  of  (5.28)  for  u2  gives  the  equation 


1  * 
7  E 

z  i,k=  1 


d2u , 


Qik 


3x;  dTf. 


/  (*) 


(5.30) 


And  this  is  the  diffusion  equation  which  defines  the  limiting  behavior  of  the  system  (5.11)  in  the 
macroscopic  £ -scale  in  the  limit  as  e  — ►  0. 

We  can  justify  the  asymptotic  analysis  by  using  energy  estimates  or  probabilistic  methods 
as  in  (Bensoussan,  Lions,  and  Papanicolaou  1978).  (See  also  Kunnemann  1983).)  We  shall  omit 
this  analysis  here. 


5.3.  Summary- 


Returning  to  the  original  problem  (5.4)  for  the  evolution  of  thermal  energy  on  a  microscopic 
scale,  we  have  shown  that  the  thermal  density  u((f  ,x  )  — *•  u0(t  ,x)  as  e  — > -  0  (in  an  appropriate 
norm)  where 


du0  i  d  du0 

~dT  =  7  .  p  ,  9ii  dx;  dx , 

*  ,  J  =  1  *  J 


with 

x  i 

9,7  =  -  ~r  Yj  {V,+  [  ai(y)Xk(y)}  +  Vk  [  a>(y)  Xi(y)  ] } 

*  k 

with  the  correctors  Xk  >  ^  —  1,2, ...,d  ,  given  by 
i 

£  vr  [  ady)  v,-+  x.-(i/)  1  =  -  [  ak(y  -  ek )  -  «k(.y)} 

*  =1 

k  —  1,2,  •  ■  •  ,  d 


(5.31) 


(5.32) 


(5.33) 
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To  compute  the  limiting  “homogenized”  model  (5.31),  one  must  solve  the  system  (5.33)  (numeri¬ 
cally)  and  then  evaluate  the  average  (5.32). 

The  fact  that  the  original  problem  (5.4)  is  “parabolic”  (i.e.,  it  describes  a  jump  random  pro¬ 
cess),  enables  us  to  exploit  the  associated  probabilistic  structure  to  anticipate  and  structure  the 
analysis.  In  this  way  we  can  anticipate  that  the  limit  problem  will  involve  a  diffusion  process.  In 
fact,  the  arguments  used  are  entirely  analytical0  and  the  limiting  diffusion  (5.31)  is  constructed  in 
a  systematic  way.  It  is  not  postulated. 


Appendix:  An  operational  calculus  for  homogenization 

The  techniques  used  to  treat  the  one  dimensional  example  are  based  on  an  operational  cal¬ 
culus  for  homogenization  developed  in  (Bensoussan,  Lions,  and  Papanicolaou  1978).  It  is 
worthwhile  to  outline  the  essential  features  of  this  operational  system. 

Consider  the  problem  (using  the  convention  that  repeated  indices  are  summed) 


A<  uc  =  - 


3 


Ky(-) 


x ,  du ' 


dxi  '  ''  c  Uxj  '  e 

X  6  O  ,  U((X  )  |  d0  =  <t> 


]  +  a  „(— )  «'=/(*) 


(1) 


where  O  is  a  smooth  domain  in  IR"  ,  30  is  its  boundary;  and  we  assume 


[  a<jiy) )  >  0  (positive  definite  matrix) 

a0  >  «  >  0  (2) 

a,y  (  .  )  ,  a  0(  •  )  periodic  with  period  1 . 

We  are  interested  in  approximating  u  c  in  the  limit  as  f  -+  0+. 

As  in  the  example,  we  introduce  y  —  x  /c  and  the  ansatz 

uc(x)  —  uc{x  )  =  u  0(x  ,—)  +  (  «  i(x  ,—)  +  ■  ■  ■  (3) 

e  e  e 

Probabilistic  arguments  can  be  used  (Bensoussan,  Lions,  and  Papanicolaou  1978,  Chapter  3);  and  they  have  some 
advantages. 
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Making  the  change  of  coordinates 


d  d  1_  _d_  _  *i_ 

dx{  dx;  e  dy{  '  e 

we  have 


A  (  =  ——  A  !  +  —  A  2  +  A  3 

e2  e 


with 


d  .  .  .  d 

1  aij(y ) 


dy, 


d  d 

a;j(y ) 


dy, 

d 


dyi 


dxj 


drj 


[  a<j(y)  ] 

9y}- 


a  d  r  /  \  5  1  /  s 

As=  -  -T7—  [  a;s(y  )  — —  ]  +  a0(y) 


dxj 


dij 


Substituting  (3)  and  (5)  into  (l)  and  equating  coefficients  of  like  powers  of  e,  we  have 


(e  2)  A  ,  u  o  =  0 
(e"‘)  A  ,  u  j  +  A2u0  =  0 
(e)  A|  #2  +  A  %  v  i  ~h  A  3  u  0  =  f 
(«‘)  +  A  2  u  2  +  A3ul  =  0 


(4) 


(4) 


(5) 


(6) 


Solvability  of  this  system  depends  on  the  operator  A}.  The  assumptions  in  (2)  mean  that 
the  system 

A^{y)  =  g  (y),  y  €  Y  (7) 

where  Y  is  the  unit  torus  in  ]Rn  and  g  (  .  )  is  periodic,  has  a  solution  unique  up  to  an  additive 
constant  if  and  only  if 

f  g  (y  )  dy  ==  0  (g) 

r 

Assumption  (2)  means  the  null  space  of  A,  is  spanned  by  the  constant  function;  and  (8)  means 
the  data  g  ( .  )  is  orthogonal  to  the  null  space  of  A  i.e.,  a  Fredholm  Alternative  holds  for  A  j. 
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Hence,  A  x  u0  =  0  implies 


«o(z>y )  =  u0(x). 


0=J41«1  +  A2  u0 

d  .  du^x.y)  d  du0(x) 

7—  [  «.•/  (y  )  a- ]  -  -a—  [  (y  )  —5- 


may  be  solved  by  taking 


u  \(x  ,y )  —  -  x3(y)  -  u0(x)  +  «i(*) 

azy 


where  the  corrector  x*  (y)  satisfies 


d  ,  ,  ,  d  { 

7: —  i  «,-y(y)  a - X  (2/) 

dy- 


The  right  side  of  (12)  satisfies  the  centering  condition  (8);  and  so,  (12)  has  a  periodic  solution, 
unique  up  to  an  additive  constant. 

Continuing  in  the  same  fashion,  the  third  equation  in  (6)  has  a  solution  «2  if 

f  { -A  2  u  1  +  A3  u  0  -  f  }  dy  =  0.  (13) 

Y 

Substituting  for  u1  from  (ll)  (12)  and  A3  from  (5),  the  solvability  condition  (13)  reduces  to 

d2u  0(x  ) 

-  Qij  -X — a -  +  Ml  ao)  uo (x)=  f  {x)  (14) 

ox i  ox  j 

uo(x)  I  r  =  <t> 


where 


<Iij  =  M[  a{j-  ]  -  M[  a,k  }  .  (15) 

The  role  of  the  corrector  in  the  approximation  is  clear  from  this.  A  simple  argument  similar  to 
the  one  used  for  the  one  dimensional  case  in  section  2  shows  that  the  matrix  {  q ,-y  }  is  non¬ 


negative  definite. 
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